% $Id: parallel1.tex 204 2009-04-20 01:59:36Z jsibert $
%
% Author: David Fournier
% Copyright (c) 2008 Regents of the University of California
%

\def\kerntwo{\kern -4pt}

\mysection{Getting Started}
For the purposes of this section a cluster of computers is
simple a collection of computers joined by a network. They may be
something exotic like a Beowolf cluster with sophistcated
network hardware or just a simple local network of machines
joined by one hub. As long as the amount of data that needs to be passed
from one computer to another is small one can successfully carry out
parallel computing even over a simple slow network.

At present the parallel distibuted computing module for \ADM\
uses the PVM (parallel virtual machine) software. 
To use the module you must install the PVM software on each computer 
in your cluster. There are installshield versions for MS operating systems
as well as various formats for Linux and other Unix systems.

Assuming that the software is installed on the various computers in your cluster
you must have an account on each system. ADM\ parallel processing
uses the master/slave model where one computer runs the master program and
there are one or more slave programs running on the computers in the
cluster. Note that more than one slave process can be run on each computer
and slave processes can be run on the same computer that is running the
master process.

\mysection{Using different compilers with PVM on Win32}
The PVM package on Win32 is compiled using Visual C++ and
the supplied libraries are Visual C++ libraries. While
many people have reported success compiling the PVM
source using other compilers we decided that the easiest
way to support all the compilers used with \ADM\ would
be to produce a single DLL that could be used by other
compilers on WIN32. A DLL is a library which is linked to
the users application at run time. To use a DLL it is necessary to link
an associated library, the import library, at link time.
A technical difficulty is that Borland and MS C++ use different
formats for their import libraries. The import libary for Borland
is named borpvm3.lib and is included with the other libraries for
\ADM\ in the various Borland binaries. The DLL is contained
in the bin subdirectory of the distribution and is named libpvm3.dll.

\mysection{A simple example}
Assume that we have $n$ $m\times m$ matrices $M_i$ indexed by $i$.
Let the $m\times m$  matrix $S$ be given by
$$S=\sum_{i=1}^n (M_i+\lambda_i \hbox{\rm I})^{-1}$$
where I denoetes the $m\times m$ identity matrix.
The problem is to minimize the squared norm of $S$, that is the sum of the
squares of the elements of $S$ with respect to the $n$ parameters~$\lambda_i$.

The (unparallelized) code for the example is
\beginexample

DATA_SECTION
  int n
  int m
 !! m=10;
 !! n=100;
 !! random_number_generator rng(131);
  3darray C(1,n,1,m,1,m)
 !! for (int i=1;i<=n;i++) C(i).fill_randn(rng);
PARAMETER_SECTION
  init_vector lambda(1,m)
  init_number kludge(2)
  3darray M(1,n,1,m,1,m)
  3darray N(1,n,1,m,1,m)
  matrix S(1,m,1,m);
  objective_function_value f
PROCEDURE_SECTION
   int i;
   for (i=1;i<=n;i++)
   {
     M(i)=C(i);
     for (int j=1;j<=m;j++)
     {
       M(i)(j,j)+=lambda(i);
     }
   }
   // compute the inverse of each matrix. 
   // this is the code which will be parallelized
   for (i=1;i<=n;i++)
   {
     N(i)=inv(N(i));
   }

   S.initialize();
   for (i=1;i<=n;i++)
   {
     S+=N(i);
   }
   f=norm2(N);
   f+=.0001*norm2(lambda);
   f+=square(kludge);
  
TOP_OF_MAIN_SECTION
  
  long int arrmblsize=5000000;
\endexample
This example shows the simplest kind of code to parallelize.
There is an outer loop and the computations for each pass
through the loop are independent of each other.
\beginexample
   for (i=1;i<=n;i++)
   {
     N(i)=inv(N(i));
   }
\endexample
Here is the code for the master program for the parallelized computations.
\beginexample
DATA_SECTION
  int n
  int m
  ivector n_per_host(1,ad_nhost)
  ivector min_n(1,ad_nhost)
  ivector max_n(1,ad_nhost)
 !! m=5;
 !! n=10;
 !! random_number_generator rng(131);
  3darray C(1,n,1,m,1,m)
 LOC_CALCS
  int i;
  for (i=1;i<=n;i++) C(i).fill_randn(rng);
  n_per_host=n/ad_nhost;
  int remainder=n%ad_nhost;
  for (i=1;i<=remainder;i++) n_per_host(i)++;
  min_n(1)=1;
  max_n(1)=n_per_host(1);
  for (i=2;i<=ad_nhost;i++) 
  {
    min_n(i)=min_n(i-1)+n_per_host(i-1);
    max_n(i)=min_n(i)+n_per_host(i)-1;
  }
  // *********  begin master constant send block  *************
  for (i=ad_nhost;i>=1;i--) 
  {
    int bufid = adpvm_master_cinitsend( PvmDataRaw );
    pvm_pack(m); // send m to slaves
    pvm_pack(C.sub(min_n(i),max_n(i))); // send C to slaves
    adpvm_master_csend(ad_stid[i-1]);
  }
  // *********  end constant send block  *************
 END_CALCS
PARAMETER_SECTION
  init_vector x(1,n)
  init_number kludge(2)
  vector f0(1,ad_nhost)
  3darray M(1,n,1,m,1,m)
  3darray sumN(1,ad_nhost,1,m,1,m)
  matrix S(1,m,1,m);
  objective_function_value f
PROCEDURE_SECTION
   int i;
   send_data(); // send x to slaves
   receive_data(); // get sumN back
   
   S.initialize();
   for (i=1;i<=ad_nhost;i++)
   {
     S+=sumN(i);
   }
   f=norm2(S);
   f+=.001*norm2(x);
   f+=kludge*kludge;
  
FUNCTION send_data
  int i;
  for (i=ad_nhost;i>=1;i--) 
  {
  // *********  begin master variable send block  *************
    int bufid = adpvm_master_vinitsend( PvmDataRaw );
    adpvm_pack(x(min_n(i),max_n(i))); // send x to slaves
    adpvm_master_vsend(ad_stid[i-1]);
  // *********  end variable send block  *************
  }

FUNCTION receive_data
  int i;
  
  double dd=0.0;
  for (i=ad_nhost;i>=1;i--) 
  {
  // *********  begin master variable receive block  *************
    adpvm_master_vrecv(ad_stid[i-1]);  // get the values from slave
    adpvm_unpack(sumN(i));
    adpvm_master_end_vreceive();  // have got all the values
  // *********  end variable receive block  *************
  }
 
GLOBALS_SECTION

  #include "adpvm2.h"

  struct pvmhostinfo *ad_hostp=NULL;
  int ad_nhost=0;
  ivector ad_stid;

TOP_OF_MAIN_SECTION
  
  long int arrmblsize=5000000;
  cout << "calling load library" << endl;
  if (load_adpvm_library() < 0)
  {
    cerr << "error loading pvm library" << endl;
    exit(1);
  }
  int result, check, i, narch;
  pvm_setopt(PvmRoute, PvmRouteDirect);  /* channel for communication */
  /* get and display configuration of the parallel machine */
  pvm_config( &ad_nhost, &narch, &ad_hostp );  /* get configuration */
  printf("I found the following hosts in your virtual machine\n");
  for (i = 0; i < ad_nhost; i++)
  { 
     printf("    %s\n", ad_hostp[i].hi_name);
  }
  ad_stid.allocate(0,ad_nhost);
  
  for (i=0; i<ad_nhost; i++)				/* spawn processes on */			
  {						/* all physical machines */
    check=pvm_spawn("pslave1", 0,PvmTaskHost /* | PvmTaskDebug */ ,ad_hostp[i].hi_name, 1, &(ad_stid(i)));
    if (!check) 
       printf("Couldn't start process on %s\n",ad_hostp[i].hi_name); 
    else
       printf("Started process on %\n",ad_hostp[i].hi_name); 
  }
FINAL_SECTION
  int i;
  // clean up slaves -- this stops all the slave processes
  for (i=0; i<ad_nhost; i++) pvm_kill(ad_stid[i]);
\endexample
The code for sending data to the slave processes and getting it back has been 
put into two functions. Here is the code to send the data. At this point
{\tt ad\_nhost} is the number of machines in the cluster and there is one slave
process running on each member of the cluster.

\beginexample
FUNCTION send_data
  int i;
  for (i=ad_nhost;i>=1;i--) 
  {
  // *********  begin master variable send block  *************
    int bufid = adpvm_\setcolor\cmykRed\kerntwo master_v\setcolor\cmykBlack\kerntwo initsend(\setcolor\cmykRed PvmDataRaw\setcolor\cmykBlack); }
   \setcolor\cmykRed 
    adpvm_pack(x(min_n(i),max_n(i))); // send x to slaves
   \setcolor\cmykBlack
    adpvm_\setcolor\cmykRed\kerntwo master_v\setcolor\cmykBlack\kerntwo send(ad_stid[i-1]);
  // *********  end variable send block  *************
  }
\endexample
All the element of the send block which are subject to change are
printed in red. If the program is a slave program 
   \setcolor\cmykRed master \setcolor\cmykBlack
 is changed to
   \setcolor\cmykRed slave \setcolor\cmykBlack\kerntwo.
If the send block is a constant send block 
   \setcolor\cmykRed v \setcolor\cmykBlack
is changed to 
   \setcolor\cmykRed c \setcolor\cmykBlack\kerntwo.
To begin the process of sending data one first calls the  {\tt adpvm\_master\_vinitsend}
function. This initializes a data buffer to hold the data prior
to sending them. The {\tt PvmDataRaw} option means that the data are not tranlated
into a format reconizeable my computers that use incompatible data
formats but are sent ``as is''. This is the faster option.
The line {\tt pvm\_pack(M.sub(min\_n(i),max\_n(i)));} put the {\tt sub dvar3\_array}
which is determined by the integers {\tt min\_(i)} and {\tt max\_(i)} into the buffer.
the {\tt sub} member function has been defined to produce sub objects of the
\ADM\ container objects. This is very convenient for communication in
parallelized code.
Now we are ready to send the data. This is accomplished with the
    {\tt adpvm\_master\_vsend((*ad\_stid)[i-1]);} 
  function. Here {\tt ad\_stid} is a
pointer to an {\tt ivector} which contains the identifier for
the i'th slave process in its i-1'th component.

The code for the slave process is not an \ADM\ program. It is a C++ progam which
uses the Autodif libraries which underly \ADM\ as well as the ADPVM module.
\beginexample

#include "adpvm2.h"

main(int argc,char * argv[])
{
  gradient_structure gs(25000000);
  time_t now;
  char name[12], buf[60];
  int ptid;

  load_adpvm_library();
  pvm_setopt(PvmRoute, PvmRouteDirect);
  
  d3_array C;
  int m;
  // ***************  begin slave constant receive block *************************
   ptid=pvm_parent();
   adpvm_slave_crecv(ptid);
   pvm_unpack(m);  
   pvm_unpack(C); 
   adpvm_slave_end_creceive();
  // ***************  end receive block ***********************************
  reset_gradient_stack();
  do
  {
    dvar_vector x;
    gradient_structure::no_derivatives=0;

    // ***************  begin slave variable receive block ***********************************
     ptid=pvm_parent();
     adpvm_slave_vrecv(ptid);
     adpvm_unpack(x);  // see if derivatives are
     adpvm_slave_end_vreceive();
    // ***************  end receive block ***********************************
  
     int nmin=C.indexmin();
     int nmax=C.indexmax();
     dvar_matrix M(1,m,1,m);
     dvar_matrix N(1,m,1,m);
     N.initialize();

     dvariable f0=0.0;
     for (int i=nmin;i<=nmax;i++)
     {
       M=C(i);
       for (int j=1;j<=m;j++)
       {
         M(j,j)+=x(i);
       }
       N+=inv(M);
     }
  
    // ***************  begin slave variable send block ***********************************
    ptid=pvm_parent();
    cout << " starting send " << endl;
    adpvm_slave_vinitsend(PvmDataDefault);/* allocate message buffer */	
    adpvm_pack(N);
    adpvm_slave_vsend(ptid);/* send buffer to master */
    cout << " finished send " << endl;
    // ***************  end send block ***********************************

    slave_gradcalc();
  }
  while(1);
}
\endexample

The code to receive the constant data from the parent program is
\beginexample
   ptid=pvm_parent();
   adpvm_slave_crecv(ptid);
   pvm_unpack(m);  
   pvm_unpack(C); 
   adpvm_slave_end_creceive();
\endexample

The code to receive the variable data from the parent program is
\beginexample
   ptid=pvm_parent();
   adpvm_slave_vrecv(ptid);
   adpvm_unpack(x);  // see if derivatives are
   adpvm_slave_end_vreceive();
\endexample
The first line {\tt int ptid=ppvm\_parent();} get a handle to the parent process
so that the slave process knows where to get its data from. The line
     { \tt adpvm\_slave\_vrecv(ptid);} receives the data and puts it into a buffer.
The line
    {\tt  pvm\_unpack(gradient\_structure::no\_derivatives);} put the value of the 
the flag into {\tt  gradient\_structure::no\_derivatives);} which turns
derivative calculations off or on as desired. 
The line {\tt  pvm\_unpack(N);}
puts the values of the {\tt dvar3\_array} into the {\tt dvar3\_array N}.
If {\tt N} is unallocated this function will allocate the right sizes for 
accepting the data that are sent. This means that one does not need to
know or calculate in advance the shape of the container object. One can simple
declare a {\tt dvar3\_array N} and everything works.

The process of sending data to the slaves has been broken up into two parts. 
Once the
{\tt d3\_array C} has been calculated a send block is set up and the parts of
{\tt C} are shipped off to the slave processes. This speeds up the 
calculations over a slow network a lot since  the larger data structure
only needs to be sent once.
\mysection {Understanding the control of the flow in a
parallelized program}

Each slave process is responsible for calculating the derivative
contributions  which correspond to the calculations it has
carried out. However all of this must occur in the opposite
order of the original calculations.

\begin{displaymath}
\xymatrix @+.85in{
{\hbox{\rm START}} \ar[d] \\
{\hbox{\rm A}}\ar^{\vbox{\hbox{\rm OTHER}\hbox{\rm CALCS}}}[d]  \ar[r]^{\hbox{\rm TO SLAVE}} & {\hbox{\rm B}} \ar[d] \\
{\hbox{\rm C}}\ar[d] & {\hbox{\rm D}} \ar[l]^{\hbox{\rm FROM SLAVE}} \ar@/_.3in/_{\vbox{\hbox{\rm LOOP}\hbox{\rm BACK}}}[u]<-.05in>\\
{\hbox{\rm END}\ar@/^.5in/^{\vbox{\hbox{\rm LOOP}\hbox{\rm BACK}}}[uuu]<.1in>} 
}
\end{displaymath}
\ADM\ sets up the derivative calculations for you automatically but
it needs a few hints. For example when the slave process receives
its data at B it does not know whether derivatives are being calculated.
If they are not being calculated the slave process should not send
the derivative information back to the master process when it
calculates the derivatives.  To let the slave process know whether 
derivatives are being calculated the first thing passed from the
master to slave process should be the integer flag named
{\tt gradient\_structure::no\_derivatives}. If this flag is non zero then
derivative calculations are turned off.
Note that the value of this flag is passed as the first thing in the
parameters passed in the \PS. It will not work to
pass it with the constant objects passed in the \DS.
%}
%\end{display}


% \beginexample
%   ad_constant_flag=1;
%   for (i=ad_nhost;i>=1;i--) 
%   {
%     int bufid = pvm_initsend( PvmDataRaw );
%     pvm_pack(m); // send m to slaves
%     pvm_pack(C.sub(min_n(i),max_n(i))); // send C to slaves
%     adpvm_send((*ad_stid)[i-1],1);
%   }
%   ad_constant_flag=0;
% \endexample
% The code for the slave processes is
% \beginexample
% #include <stdio.h>
% extern "C" {
% #include <pvm3.h>
% }
% #include <admodel.h>
% #include <adpvm.h>
% #if !defined(linux)
% #  include <windows.h>
% #endif
% #define  __declspec(dllexport) 
% 
% extern int ad_constant_flag;
% 
% main(int argc,char * argv[])
% {
%   gradient_structure gs(25000000);
%   time_t now;
%   char name[12], buf[60];
%   int ptid;
% 
%   load_adpvm_library();
%   pvm_setopt(PvmRoute, PvmRouteDirect);
%   
%   d3_array C;
%   int m;
%   // ***************  begin receive block ***********************************
%    ad_constant_flag=1;
%    ptid=pvm_parent();
%    adpvm_recv(ptid,-1);
%    pvm_unpack(m);  // see if derivatives are
%    cout << "m = " << m << endl;
%    pvm_unpack(C); 
%      adpvm_end_receive();
%    ad_constant_flag=0;
%   // ***************  end receive block ***********************************
%   do
%   {
%     dvar_vector x;
% 
%     // ***************  begin receive block ***********************************
%    ptid=pvm_parent();
%      adpvm_recv(ptid,-1);
%      pvm_unpack(gradient_structure::no_derivatives);  // see if derivatives are
%      pvm_unpack(x);  // see if derivatives are
%      adpvm_end_receive();
%     // ***************  end receive block ***********************************
%   
%      int nmin=C.indexmin();
%      int nmax=C.indexmax();
%      dvar_matrix M(1,m,1,m);
%      dvar_matrix N(1,m,1,m);
%      N.initialize();
% 
%      dvariable f0=0.0;
%      for (int i=nmin;i<=nmax;i++)
%      {
%        M=C(i);
%        for (int j=1;j<=m;j++)
%        {
%          M(j,j)+=x(i);
%        }
%        N+=inv(M);
%      }
%   
%     // ***************  begin send block ***********************************
%    ptid=pvm_parent();
%     cout << " starting send " << endl;
%     pvm_initsend(PvmDataDefault);/* allocate message buffer */	
%     pvm_pack(N);
%     adpvm_send(ptid, 2);/* send buffer to master */
%     cout << " finished send " << endl;
%     // ***************  end send block ***********************************
% 
%     cout << "starting slave gradcalc " << endl;
%     slave_gradcalc();
%     cout << "finished slave gradcalc " << endl;
%   }
%   while(1);
% }
% 
% \endexample
% The code for the extra receive block is below. As in the send block for the
% master process.  The {\tt ad\_constant\_flag} must be set so that the proper code is
% produced.
% \beginexample
%    ad_constant_flag=1;
%    ptid=pvm_parent();
%    adpvm_recv(ptid,-1);
%    pvm_unpack(m); 
%    pvm_unpack(C); 
%    adpvm_end_receive();
%    ad_constant_flag=0;
% \endexample
% With these modifications the parallelized code ran about twice as fast
% the unparallelized code.
